# For use in conjunction with S4-classes based modules in cRy.
# The word "dot" has been adopted because all functions here start with a ".", so
# to be invisible.


.dcl <- function(a,b,c,aa,bb,gg)
{
# Input: cell parameters. Output: values useful to all crystallographic
# calculations. These are: 
# 1) sa = sin(alpha), sb = sin(beta), sc = sin(gamma)
# 2) ca = cos(alpha), cb = cos(beta). cc = cos(gamma)
# 3) sides of reciprocal cell: ar = a*, br = b*, cr = c*
# 4) sines of angles of reciprocal cell: sar = sin(alpha*), sbr = sin(beta*), scr = sin(gamma*)
# 5) cosines of angles of reciprocal cell: car = cos(alpha*), cbr = cos(beta*), ccr = cos(gamma*)
# 6) Volume of unit cell: V
 aa <- aa*pi/180
 bb <- bb*pi/180
 gg <- gg*pi/180
 sa <- sin(aa)
 sb <- sin(bb)
 sc <- sin(gg)
 ca <- cos(aa)
 cb <- cos(bb)
 cc <- cos(gg)

 # To avoid NaN generated by rounding off errors, use for cell-derived quantities formulas
 # derived previously by computationa crystallographers
 sang <- 0.5*(aa+bb+gg)
 V2 <- sqrt(sin(sang-aa)*sin(sang-bb)*sin(sang-gg)*sin(sang))
 V <- 2*a*b*c*V2
 ar <- b*c*sa/V
 br <- a*c*sb/V
 cr <- a*b*sc/V
 car <- (cb*cc-ca)/(sb*sc)
 cbr <- (ca*cc-cb)/(sa*sc)
 ccr <- (ca*cb-cc)/(sa*sb)
 sar <- sqrt(1-car*car)
 sbr <- sqrt(1-cbr*cbr)
 scr <- sqrt(1-ccr*ccr)
 l <- c(sa,sb,sc,ca,cb,cc,ar,br,cr,sar,sbr,scr,car,cbr,ccr,V)
 names(l) <- c("SIN_ALPHA","SIN_BETA","SIN_GAMMA","COS_ALPHA","COS_BETA","COS_GAMMA","A*","B*","C*","SIN_ALPHA*","SIN_BETA*","SIN_GAMMA*",
               "COS_ALPHA*","COS_BETA*","COS_GAMMA*","V")

 return(l)
}

#
.crystal_system <- function(gn)
{
 # Given the space group number (gn), returns the crystal system
 # bs = 1 TRICLINIC
 # bs = 2 MONOCLINIC
 # bs = 3 ORTHOROMBIC
 # bs = 4 TETRAGONAL
 # bs = 5 CUBIC
 # bs = 6 HEXAGONAL
 # bs = 7 TRIGONAL
 if (gn >=   1 & gn <=   2) bs <- 1
 if (gn >=   3 & gn <=  15) bs <- 2
 if (gn >=  16 & gn <=  74) bs <- 3
 if (gn >=  75 & gn <= 142) bs <- 4
 if (gn >= 143 & gn <= 167) bs <- 7
 if (gn >= 168 & gn <= 194) bs <- 6
 if (gn >= 195 & gn <= 230) bs <- 5
 bs_name <- c("TRICLINIC","MONOCLINIC","ORTHOROMBIC","TETRAGONAL","CUBIC","HEXAGONAL","TRIGONAL")

 return(bs_name[bs])
}

# Given cell parameters, return the inverse of cell orthogonalization matrix (First choice in Giacovazzo's book)
.triclinic_to_orthogonal_01 <- function(a,b,c,aa,bb,cc)
{
 lp <- .dcl(a,b,c,aa,bb,cc)
 m_1 <- matrix(c(a,0,0,b*lp[6],b*lp[3],0,c*lp[5],-c*lp[2]*lp[13],1/lp[9]),nrow=3,ncol=3,byrow=TRUE)

 return(m_1)
}

# Given cell parameters, return the inverse of cell orthogonalization matrix (MOSFLM choice, second choice in Giacovazzo's book)
.triclinic_to_orthogonal_02 <- function(a,b,c,aa,bb,cc)
{
 lp <- .dcl(a,b,c,aa,bb,cc)
 m_1 <- matrix(c(1/lp[7],-lp[15]/(lp[7]*lp[12]),a*lp[5],0,1/(lp[8]*lp[12]),b*lp[4],0,0,c),nrow=3,ncol=3,byrow=TRUE)

 return(m_1)
}

#
# Transform fractional into orthonormal coordinates
.frac_to_orth <- function(xyzf,a,b,c,aa,bb,cc,ochoice=1)
{
 # a,b,c = cell sides in angstroms; aa,bb,cc = cell angles in degrees
 # The parameter ochoice controls which convention is being used to
 # collocate the cell in an orthonormal cartesian frame.
 # ochoice = 1: X axis along a; Y axis normal to a, in the (a,b) plane;
 #              Z axis normal to X and Y (and therefore parallel to
 #              c*).
 # ochoice = 2: this is also called "Cambridge setting". The X axis is
 #              along a*; the Y axis lies in the (a*,b*) plane; the Z
 #              axis is, consequently, along c. 
 # xyzf is a vector, matrix or data frame of fractional crystal coordinates.
 # If matrix or data frame it needs to have 3 columns. 
 # This function returns a data frame with 3 columns, the cartesian orthonormal coordinates.

 # Check xyzf object is either vector, matrix or data.frame
 if (!is.vector(xyzf) & !is.matrix(xyzf) & !is.data.frame(xyzf)) stop("Input object is not of a valid type (vector, matrix or data.frame)")

 # If xyzf is a vector, turn it into a matrix
 if (is.vector(xyzf)) dim(xyzf) <- c(1,length(xyzf))

 # If xyzf has less or more than 3 columns stop
 if (dim(xyzf)[2] != 3) stop("Input object has more or less than 3 columns")

 # If xyzf is a data.frame turn it into a matrix
 if (is.data.frame(xyzf)) xyzf <- as.matrix(xyzf)

 # Orthogonalization matrix
 if (ochoice == 1) M_1 <- .triclinic_to_orthogonal_01(a,b,c,aa,bb,cc)
 if (ochoice == 2) M_1 <- .triclinic_to_orthogonal_02(a,b,c,aa,bb,cc)

 # Transformed coordinates
 #x <- av[1]*xf+bv[1]*yf+cv[1]*zf
 #y <- av[2]*xf+bv[2]*yf+cv[2]*zf
 #z <- av[3]*xf+bv[3]*yf+cv[3]*zf
 xyz <- xyzf%*%M_1 

 # Turn matrix back into a data frame
 xyz <- as.data.frame(xyz)
 colnames(xyz) <- c("x","y","z")

 return(xyz)
}

#
# Transforms orthogonal into fractional coordinates
.orth_to_frac <- function(xyz,a,b,c,aa,bb,cc,ochoice=1)
{
 # Does the inverse job of "frac_to_orth"

 # Check xyz object is either vector, matrix or data.frame
 if (!is.vector(xyz) & !is.matrix(xyz) & !is.data.frame(xyz)) stop("Input object is not of a valid type (vector, matrix or data.frame)")

 # If xyz is a vector, turn it into a matrix
 if (is.vector(xyz)) dim(xyz) <- c(1,length(xyz))

 # If xyz has less or more than 3 columns stop
 if (dim(xyz)[2] != 3) stop("Input object has more or less than 3 columns")

 # If xyz is a data.frame turn it into a matrix
 if (is.data.frame(xyz)) xyz <- as.matrix(xyz)

  # Orthogonalization matrix
 if (ochoice == 1) M_1 <- .triclinic_to_orthogonal_01(a,b,c,aa,bb,cc)
 if (ochoice == 2) M_1 <- .triclinic_to_orthogonal_02(a,b,c,aa,bb,cc)

 # Inverse matrix
 M <- solve(M_1)

 # Transform coordinates
 xyzf <- xyz%*%M

 # Turn matrix back into a data frame
 xyzf <- as.data.frame(xyzf)
 colnames(xyzf) <- c("xf","yf","zf")

 return(xyzf)
}

.d_hkl <- function(h,k,l,a,b,c,aa,bb,cc)
{
 # Given Miller indices and cell parameters, this function returns
 # resolution corresponding to the specific Miller indices.

 aa <- aa*pi/180
 bb <- bb*pi/180
 cc <- cc*pi/180
 top <- 1-(cos(aa))^2-(cos(bb))^2-(cos(cc))^2+2*cos(aa)*cos(bb)*cos(cc)
 b1 <- h^2*(sin(aa))^2/a^2
 b2 <- k^2*(sin(bb))^2/b^2
 b3 <- l^2*(sin(cc))^2/c^2
 b4 <- 2*h*k*(cos(aa)*cos(bb)-cos(cc))/(a*b)
 b5 <- 2*h*l*(cos(aa)*cos(cc)-cos(bb))/(a*c)
 b6 <- 2*k*l*(cos(bb)*cos(cc)-cos(aa))/(b*c)
 d2 <- top/(b1+b2+b3+b4+b5+b6)
 return(sqrt(d2))
}
